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ABSTRACT 


Context. Atmospheric properties of exoplanets can he constrained with transit spectroscopy. At low spectral resolution, this technique 
is limited by the presence of clouds. The signature of atomic sodium (Nai), known to he present above the clouds, is a powerful probe 
of the upper atmosphere, where it can be best detected and characterized at high spectral resolution. 

Aims. Our goal is to obtain a high-resolution transit spectrum of HD 189733b in the region around the resonance doublet of Nai at 
589 nm, to characterize the absorption signature previously detected from space at low resolution. 

Methods. We analyze archival transit data of HD 189733b obtained with the HARPS spectrograph (K = 115 000) at ESO 3.6-meter 
telescope. We perform differential spectroscopy to retrieve the transit spectrum and light curve of the planet, implementing corrections 
for telluric contamination and planetary orbital motion. We compare our results to synthetic transit spectra calculated from isothermal 
models of the planetary atmosphere. 

Results. We spectrally resolve the Nai D doublet and measure line contrasts of 0.64±0.07% (D2) and 0.40±0.07% (Dl) and FWHMs 
of 0.52 ± 0.08 A. This corresponds to a detection at the 10-cr level of excess of absorption of 0.32 ± 0.03% in a passband of 2 x 0.75 A 
centered on each line. We derive temperatures of 2 600 ± 600 K and 3270 ± 330 K at altitudes of 9 800 ± 2 800 and 12 700 ± 2 600 km 
in the Nai Dl and D2 line cores, respectively. We measure a temperature gradient of ~ 0.2 K km“' in the region where the sodium 
absorption dominates over the haze-absorption, from comparison with theoretical models. We also detect a blueshift of 0.16 ± 0.04 A 
(4 (t) in the line positions. This blueshift may be due to winds blowing at 8 ± 2 km s“’ in the upper layers of the atmosphere. 
Conclusions. We demonstrate the relevance of studying exoplanet atmospheres with high-resolution spectrographs mounted on 4- 
meter-class telescopes. Our results pave the way towards in-depth characterization of physical conditions in the atmospheres of many 
exoplanetary systems with future spectrographs such as ESPRESSO on the VET or HiReS and METIS on the E-ELT. 

Key words. Planetary Systems - Planets and satellites: atmospheres, individual: HD 189733b - Techniques: spectroscopic - 
Instrumentation: spectrographs - Methods: observational 


1. Introduction 


After two decades and over a thousand exoplanet detections, we 
have entered an era of characterization of these remote worlds. 
More and more exoplanets have precise mass and radius mea¬ 
surements thanks to simultaneous radial-velocity and transit ob¬ 
servations, providing us with estimates of mean densities and 
constraints on possible bulk compositions. Meanwhile, detec¬ 
tions and studies of exoplanet atmospheres represent the only 
direct observable window we have on the physical and chemical 
properties of these planets. Characterization of exoplanets and 
comparative planetology are now progressing well by virtue of 
cutting-edge instrumentation. Among these observations, spec¬ 
tra of exoplanets are among the best products for in-depth com¬ 
prehension of atmospheric and surface conditions. 

Shortly after the first indirect exoplanet discoveries, early ob¬ 
servational campaigns encouraged by theoretical work ( |Marley 
et al.|199^ Seager & Sasselov|2000[ Brown|2001|l, attempted to 


detect reflected light from the day side ( Charbonneau et al.|1999| l 
or absorption through limb transmission (Moutou et al. |2001| l. 
Observations of HD 209458, the unique transiting system known 
at that time, with the Space Telescope Imaging Spectrograph 


* using observations with the Harps spectrograph from the ESO 3.6m 
installed at La Silla, Chile, under the allocated programmes 072.C- 
0488(E), 079.C-0828(A) and 079.C-0127(A). 


(STIS) on board the 2.4-m Hubble Space Telescope (HST), al¬ 
lowed [Charbonneau et ah] ( |2002| ) to detect atmospheric sodium 
(Nai) at 589 nm. This measurement was performed with a spec¬ 
tral resolution of 7? s 4/AT ~ 5 500 (or ~ 55 km s“'). Extra ab¬ 
sorption of 0.023+0.006% and 0.013+0.004% have been mea¬ 
sured during transits over spectral bins of 12 A and 38 A, re¬ 
spectively. The Nai doublet lines Dl (4589.5924 nm) and D2 
(4588.9951 nm) are not resolved in these bins. 


In addition to important detections of other atomic and 
molecular signatures (e.g. jVidal-Madjar et ah [2003 1 jPemingj 
et al.|201? I, the sodium signature in HD 209458b has remained 
one of the most robust examples of atmospheric characteriza- 
tion (jMadhusudhan et al.p014| [Heng & Showman||2014[ jPepej 


et al. 2014 1 . It has been confirmed by independent analysis of the 
same data set ( [Sing et ar]|2008a|b) l. In the early days, ground- 
based observations of transiting systems with high-resolution 
spectrographs mounted on 8-meter-class telescope were also 
quickly tested in the perspective of confirming detections and 
resolving transmission spectra, unfortunately without any pos¬ 
itive results (e.g. Narita et al. |2005| l. Indeed, even with high- 
enough signal, no good observational strategy or data analysis 
were found to overcome terrestrial atmospheric variation and 
telluric lines contamination that transmission spectra undergo. 
Later, careful data (re-)analysis made by Redfleld et al. ( |2008| l 
and [Snellen et al.|(|2008|l enabled to beat down systematics due 
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to Earth atmosphere and led to the first ground-based detec¬ 
tion of exoplanet atmospheric signature. Indeed the Nat dou¬ 
blet was detected for the first time in HD 189733b with the High 
Resolution Spectrograph (HRS; 7? ~ 60000) mounted on the 9- 
m Hobby-Eberly Telescope and confirmed in HD 209458b using 
the High Dispersion Spectrograph (HDS; 7? ~ 45 000) on the 
8-m Subaru telescope. Eurthermore, the analysis of the space- 
based detection of the sodium signature of these two planets by 
Sing et al. (2008b i, VidaTMadjar et al.] ( 201 1[ ), and Huitson et al. 
( 2012| l, demonstrated the potential to resolve the transmitted 
star light by the planet atmosphere. Despite these encouraging 
works, the presence of scattering hazes and clouds (Lecavelier 


Des Etangs et al.|2008a|b[ Pont et al.|2013| l in several exoplanets 
prevented the detection of major chemical constituents at low-to 
mid-resolution even from space (e.g. Kreidberg et al.|20T4{ Sing 
et al.|2015 1 . The whole potential of high-resolution spectroscopy 
with R ~ 100 000 applied to transiting and non-transiting system 
was shown when H 2 O and CO molecular bands were detected on 
account of cross-correlation on hundreds of resolved individual 
transitions ( Snellen et al.|201()l 2014[ Brogi et al.|2012[ Birkby 


et al.|201^ . This state-of-the-art method brought back attention 

to the importance of transit spectroscopy observations from the 
ground. 

Transit spectroscopy with ground-based telescopes is of¬ 
ten used to gauge the Rossiter-McLaughlin effect during an 
exoplanet transit. This is accomplished with high-precision 
velocimeters. As data are acquired during transit, we can 
think about measurements of exoplanet atmosphere properties 
with the same data. However, observations with simultaneous- 
calibrated or self-calibrated velocimeter are not equally valuable 
because of the use of an iodine cell. The absorption of the light 
by the iodine cell act as a contaminant. The iodine lines pre¬ 
vent measurements of several spectral features especially around 
the Nat lines. Thus the use of stabilized spectrographs with si¬ 
multaneous calibration makes possible to measure the Rossiter- 
McLaughlin effect and the transmission spectrum of an exo¬ 
planet at the same time. 

In this paper, we will describe new results obtained with 
the HARPS spectrograph by analyzing existing high-quality 
Rossiter-McLaughlin observations of HD 189733b. Eor the first 
time, features in an exoplanet atmosphere are resolved with a 4- 
m-class telescope. In Sect. we will present the observations, 
the data reduction, the telluric correction and the method to mea¬ 
sure exo-atmospheric signals. In Sect.[^ we will describe the re¬ 
sults obtained on the measure of the Nat D absorption excess in 
HD 189733b. In Sect. we will discuss the implication of our 
detection, eventually we will compare our results with theory. 


2. Observations, data reduction and methods 

2.1. HARPS observations 


HD 189733 ( Bouchy et al.|2005] l was observed with the HARPS 
echelle spectrograph (High-Accuracy Radial-velocity Planet 
Searcher) on the ESO 3.6 m telescope. La Silla, Chile. Data were 
retrieved from the ESO archive from programs 072.C-0488(E), 
079.C-0828(A) (PkMayor) and 079.C-0127(A) (PI: Lecavelier 
des Etangs). The data sets are described in Table [T] (nights 1,3,4 
were employed by |Triaud et al.| ( [2009| l to measure the Rossiter- 
McLaughlin effect). In total, four sequences of data cover at least 
partially the transit of the planet. Among those four sequences, 
two were obtained using low cadence (600 s exposures), the two 
others using high cadence (300 s exposures). One low-cadence 
sequence was affected by bad weather and missed the second 


half of the transit. Eor this reason, we decided not to take this 
sequence (2006-07-29) into account. Our analysis is based on 
the remaining three data sequences totalizing 99 spectra and 3 
planetary transits. Based on the transit ephemeris of |Agol et al.| 
( 2010| l (see Table |^, we identihed 46 spectra fully obtained in¬ 
transit (i.e. fully between the first and fourth contacts). The 53 
remaining spectra constitute our out-of-transit sample. 


2.2. Data reduction 

The HARPS observations were reduced with the version 3.5 of 
the HARPS Data Reduction Software. Spectra are extracted or¬ 
der by order and flat-helded with the daily calibration set. Each 
spectral order is blaze-corrected and wavelength-calibrated. All 
spectral orders from a given two-dimensional echelle spectrum 
are merged and resampled with a 0.01 A wavelength step into a 
one-dimensional spectrum. Elux is conserved during this step. A 
reduced HARPS spectrum covers the region between 380 nm 
and 690 nm with a spectral resolution of 7? ~ 115 000 or 
2.7 km s“'. All spectra are referred to the Solar System barycen- 
ter rest frame and wavelengths are given in the air. 

As we will see in Sect. |2.4[ the next step is to co-add the in¬ 
transit spectra on the one hand, and the out-of-transit spectra on 
the other hand to build the master-in and master-out spectra. The 
ratio of the master-in and the master-out yields the transmission 
spectrum. Then, any changes in the stellar line shape or posi¬ 
tion during the observations can create spurious signals in the 
transmission spectra, which can become false-positive of atmo¬ 
spheric signals (i.e., with similar amplitude). This is the case for 
changes in the line spread function (ESP) of the instrument. This 
is where the spectrograph design is important: HARPS is a hber- 
fed spectrograph stabilized in temperature and pressure, ensur¬ 
ing that changes in the point spread function (PSP) between two 
consecutive observations lead to negligible changes in the ESP 
at time scales of one night. Nonetheless, change in the line posi¬ 
tions are due to the radial-velocity variation of the star during the 
transit of the planet (~ 50 m s '). This must be removed before 
computing the transmission spectra. Actually, omitting this step 
for spectral shifts of 0.01 A or larger creates artificial signals. We 
simply correct from this effect by shifting the spectra to the null 
stellar radial velocity (i.e. the mid-transit stellar radial velocity) 
according to a Keplerian model of the targeted planetary system 
(see the model parameters in Table|^. Another source of change 
in the stellar line shape is the Rossiter-McLaughlin effect. As 
the transit occurs, the planet occults different parts of the stellar 
disk which have different intrinsic line shapes and shifts. Palse- 
positive features will then appear in the transmission spectrum 
even if the stellar spectra are well aligned. A way to overcome 
this difficulty is to register spectra uniformly during the transit, 
so that this effect will be averaged out during a full transit. 


2.3. Teiiuric correction 


High-resolution spectra recorded from the ground bear the im¬ 
print of Earth’s atmosphere. In the visible domain covered with 
HARPS, water vapor and molecular oxygen are the main con¬ 
tributors to this time-variable telluric contamination. Indeed, 
variation in the transmission of Earths atmosphere during a night 
depends on the airmass and on water column variations in the air. 
In addition to the contaminating water lines around the Na i dou¬ 
blet, we expect the observed spectra to also possess a signature 
of the telluric sodium ( |Vidal-Madjar et ah 2010 1 . Even if the 
telluric sodium undergoes seasonal variations and possibly does 
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Table 1. Observations log for HD 189733. Three transit events are analyzed in our work. 



Date 

# Spectra* 

Exp. Time [s] 

Airmass 

Seeing 

SNR2 

SNR3 

Program 

Analysis 

Night 1 

2006-09-07 

20 (9/11) 

900 to 600 

1.6-2.1 

0.7-1.0 

~ 170 

~ 32 

072.C-0488(E) 

Yes 

Night 2 

2007-07-19 

39 (18/21) 

300 

2.4-1.6 

0.6-0.8 

~ 110 

~ 20 

079.C-0828(A) 

Yes 

Night 3 

2007-08-28 

40 (19/21) 

300 

2.2-1.6 

0.7-2.0 

~ 100 

~ 18 

079.C-0127(A) 

Yes 

Night 4 

2006-07-29 

12 (5/7) 

600 

1.8-1.6 

0.9-1.2 

~ 140 

~ 27 

072.C-0488(E) 

No 


Notes, 'in parenthesis: the number of spectra taken during the transit (forming our in-transit sample) and outside the transit (forming our out- 
of-transit sample). ^The signal-to-noise ratio (SNR) per extracted pixel in the continuum near 590 nm. ^The SNR in the lines core of the Nai D 
doublet (this show the hrst difficulty of sodium detection in exoplanet atmospheres). 


Table 2. Adopted values for the orbital and physical parameters of HD 189733b. 


Parameters 

Symbol 

Values 

Units 

References 


Transit epoch (BJD) 

To 

2454279.436714 ± 0.000015 

days 

Agol et al. I 

2010 


Orbital Period 

P 

2.21857567 ± 0.00000015 

days 

Ago! et at. 

2010 


Planet/star area ratio 

(RpIRs? 

0.02391 ± 0.00007 


lorres et al. 

2008 


Transit duration 

h 

0.07527 ± 0.00037 

days 

iriaud et at. 

2009 

1 

Impact parameter 

b 

0.6631 ±0.0023 

R* 

lAgoi et at. 

2UTUr 

Orbital semi-major axis 

a 

0.0312 ±0.00037 

AU 

Iriaud et at. 

2009 

Orbital inclination 


85.710 ±0.024 

degrees 

Ago! et at. 

2010 1 

Orbital eccentricity 

e 

0 


hxed 



Longitude of periastron 

(1) 

90 

degrees 

fixed 



Stellar velocity semi-amplitude 

Ki 

200.56 ± 0.88 

km i“* 

Boisse et al. 

20091) 

Systemic velocity 

7 

-2.2765 ± 0.0017 

km i“* 

Boisse et al. 

2009 

Stellar mass 

M* 

0.823 ± 0.029 

Mo 

iriaud et at. 

20091 

Stellar Radius 

R* 

0.756 ±0.018 

Re 

lorres et al. 

20081 

Planet Mass 

Mp 

1.138 ±0.025 

Mj 

iriaud et ai. 

2009 


Planet Radius 

Rp 

1.138 ±0.027 

Rj 

i’orres et al. (2008) 


not follow the water absorption levels, we make the assumption 
that, within a night, this telluric sodium absorption behave and 
can be corrected as other telluric features ( [Snellen et al.||2008] 
Zhou & Bayliss|20T2 1. Furthermore, visual inspection of our tel¬ 
luric spectra did not reveal any obvious telluric sodium features. 

Telluric contamination produce ubiquitous features in trans¬ 
mission spectra. However, the telluric features can be removed 
by subtracting high-quality telluric spectra obtained on the same 
nights. As our observed spectra are referred to the Solar System 
barycentric rest frame, the telluric lines are shifted during the 
three different nights by -0.26 A, 0.09 A and -0.19 A compared 
to their rest frame. Methods describ ed in jVidal-Madjar et al. 
( 2010| l and Astudillo-Defru & Rojo ( 2013| permit us to build 
high-quality telluric spectra night hy night using the HARPS ob¬ 
servations themselves. Both methods consider that the variation 
of telluric lines follows linearly the airmass variation. This is a 
consequence of the usual hypothesis of a radiative transfer in 
a plane parallel atmosphere. It implies that apart from different 
constructions of telluric spectrum these methods yields equiva¬ 
lent telluric spectra. Similarly, we build for each night a refer¬ 
ence telluric spectrum T (d) such as 


T{X) ^T\A) = 


( 1 ) 


corresponding to a zenithal atmospheric transmission (airmass 
secz = 1). The zenithal optical depth is derived at each 
wavelength A by linear regression on the logarithm of the nor¬ 
malized flux log F(A) as a function of sec z ( [Astudillo-Defru &| 
Rojo|2013]). The computed telluric spectra (see Fig.[]J) are con 


sistent with the atmospheric transmission obtained by Hinkle 
[et al.[(|200^. Note that the telluric spectra are built with the out- 


of-transit sample only. This is important in order not to over¬ 
correct for the possible exoplanet atmosphere we are looking 
for. However, for night 2, we had to take into account the in¬ 
transit sample because the base line was insufficient to compute 
a high-quality telluric spectrum (during this night, there was no 
observation before transit). 

Taking advantage of the scaling relation between the telluric 
line strength and the airmass, we can correct spectra for telluric 
contamination. Corrections are made by rescaling all spectra as 
if they had been observed at the same airmass. To obtain this 
spectrum Fref at airmass sec Zref, we simply divided the ob¬ 
served spectra OseczW taken at airmass sec z by the telluric ref¬ 
erence spectrum to the power sec z - sec Zre/'- 

Fref(X) = ■ (2) 

We hx the reference airmass to the average airmass of the in¬ 
transit spectra. This ensures both a telluric correction and min¬ 
imal flux changes in the reduced spectra. If the quality and the 
stability of the observation nights are good enough which typ¬ 
ically depends on a constant water atmospheric content, this 
method is sufficient to correct down to the statistical noise level. 
However, the quality of the nights are generally not sufficient as 
second-order variations of telluric lines are present and telluric 
residual can be still visible in transmission spectra (see Fig.[T]and 
Sect. [2.4| l. Indeed these observation were not planned to observe 
exoplanetary atmospheres but Rossiter-McLaughlin effects. This 
is less demanding in terms of sky quality and requires shorter 
time baseline. Later in our analysis, we will perform a second 
telluric correction to remove possible telluric residuals. 
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2.4. Transmission spectrum 

Transit spectroscopy is a differential technique requiring the ac¬ 
quisition of spectra during and outside transit events. Ground- 
based observations out-of-transit (Oout) contain the star light ab¬ 
sorbed by the terrestrial atmosphere. Observations during tran¬ 
sit (Oin) additionally contain the exoplanet atmosphere trans¬ 
mission diluted into the stellar flux. Stacking telluric-corrected 
spectra in-transit (F,„) and out-of-transit (Fgut) allows to ob¬ 
tain master-in {f^m) and master-out spectra (Tout), respectively. 
Classical methods obtain transmission spectra (9? = 91 - 1 = 
'FinI'Fout - 1) of the exoplanet atmosphere by dividing night by 
night the master spectrum in-transit by the master spectrum out- 
of-transit ( |Redfield et al.|2008| l. We can then analyze exoplanet 
atmospheres features along wavelengths of interest. Such meth¬ 
ods do not consider changes in radial velocity of the planet. 

Here, we present a modified approach with respect to 
[Redfield et al.|p008) to correct for this effect. As the planet tran¬ 
sits its star, its radial velocity changes typically from -15 km s“^ 
to h- 15 km s"'. Therefore, for the sodium the planetary atmo¬ 
spheric absorption lines move from the blue-shifted to the red- 
shifted part of the stellar lines. In the optical, near the sodium 
doublet this shift is < 0.5 A during the full transit. We correct 
this effect in our analysis by shifting the planetary signal to the 
null radial velocity in the planet rest frame (i.e. the radial veloc¬ 
ity of the planet in the middle of the transit). Instead of taking the 
ratio of master spectra, we divide each spectrum in-transit by the 
master-out, then shift the residual with the computed planetary 
radial velocity and sum up all the individual transmission spec¬ 
tra. We normalized this sum to unity by a linear fit to the contin¬ 
uum outside the regions of interests. This allows us to compute 
the transmission spectra as 


J out 


Fin(A) 


2 Fouti-T) 'Planet RV shift 
out 


- 1 


(3) 


With this correction, we expect to see the exoplanetary signal at 
wavelengths corresponding to the systemic velocity rest frame 
(and not in the solar system rest frame). 

A transmission spectrum (with or without radial-velocity 
correction) corrected from the airmass effect (Sect. |23] l can still 
present some telluric line residuals. This is due to water col¬ 
umn variations in the air above the telescope. Then, we applied 
a second telluric correction to correct from telluric features in 
the transmission spectra. Because both the transmission spectra 
and the telluric spectra contain telluric features at corresponding 
wavelengths, these two spectra are correlated. To eliminate these 
correlations and the residuals of telluric lines in the transmission 
spectrum, a linear fit is made between all corresponding fluxes 
in wavelength of the transmission spectrum and the telluric ref¬ 
erence spectrum. The transmission spectrum is then divided by 
the fit solution. As only telluric residuals are correlated with the 
telluric spectrum, doing this iteratively eliminates telluric resid¬ 
uals in the transmission spectra without affecting the other part 
of the spectra ( Snellen et al.|[2008| 2010[ Khalafinejad 2013| l. 
Practically, one iteration is sufficient to tame atmospheric pollu¬ 
tion down to the photon-noise level. 


2.5. Binned atmospheric absorption depth 

In order to compare our results with previous detections of 
sodium in HD 189733b ( [Redfield et al.|20()8||Jensen et al.|20TT 
Huitson et al.||2012|l or in other exopla nets ^Charbonneau et al. 
2002| [Snellen et al.||2008| [Sing et aT][2008a[ [Langland-Shula 



5900.0 5901.0 5902.0 


A [A] 

Fig. 1. Example of the effect of the Earth atmosphere on exo¬ 
planetary transmission spectra. This shows the importance of a 
precise telluric correction. Upper panel; Calculated telluric spec¬ 
trum T(A) (top), transmission spectrum without any correction 
(middle), and transmission spectrum with the telluric correction 
(bottom) for night 2. The telluric lines are dominating all fea¬ 
tures in the uncorrected transmission spectrum. Lower panel; 
Transmission spectra with telluric corrections already performed 
for night 3. Once we apply the telluric correction described in 
Sect. [2.3[ (airmass effect) we get a transmission spectra with 
weak telluric lines residuals (top). Since residuals of telluric 
lines are still visible, we apply a second telluric correction (cor¬ 
rection of features due to water column variation, see Sect. [2.4[ l. 
This affords us to get a transmission spectrum corrected from the 
effect of Earth atmosphere (bottom). Note the different vertical 
scales. Some evident telluric lines are emphasized with vertical 
lines. The telluric lines are shifted from one night to the other 
because the spectra are referred to the Solar System barycenter 
rest frame and were not observed at the same barycentric Earth 
radial velocity. 


et al.|2009[[Wood et al.|20TT][Sing et al.|2()12[[Zhou & Bayliss 


2012 


Murgas et al.[[2014t [Nikolov et al.[[2014[ [Burton et al. 


20151, we calculate the relative absorption depth across vari¬ 


ous bins in wavelength. Taking into account the systemic ve¬ 
locity, we average the flux around the central passbands (C) 
centered on both lines of the sodium doublet. Here and later in 
our analysis, all the averages are weighted by the estimated er¬ 
rors on the reduced spectra, which are taken to be random pho¬ 
ton noise obeying Poisson statistics. We compare the integrated 
flux in the transmission spectrum 9i (C) to bins of similar band 
widths taken in the transmission spectrum continuum. Snellen 
et al. ( [2008[ ) chose adjacent control passbands on the blue (B) 
and red (R) side of the central passband. As at high-resolution 
exoplanet atmospheric lines are likely to be resolved, we pre¬ 
fer to choose absolute reference passbands outside the sodium 
doublet, but still on the two side (B and R) of the transition 
(see e.g. Charbonneau et al.|2002[ l. Relative depths due to exo- 
atmospheric absorption are then obtained by the difference of 
fluxes between the central and the reference passband. 


6(AA) = 91'(C) - 


91'(B) H-91'(R) 
2 


(4) 


As small passbands encompass only one line, we average the 
absorption depth of the two Na i D lines. When the central pass- 
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band includes the two lines, we adjusted it on the center of the 
doublet. 


2.6. Transmission iight curve 

This method consists in directly deriving the relative transmis¬ 
sion light curve as a function of time ( Charbonneau et al.|2002 
Snellen et al. 20081. Thus, the absorption excess due to the 


exoplanet atmospheric limb can be seen as a relative flux de¬ 
crease during the transit. It differs from the method presented 
in Sect. |2.4| The spectrophotometry is performed on individual 
spectra for a given spectral bin. We obtain a time-dependent in¬ 
formation but the wavelength dependence is lost. 

For each telluric-corrected stellar spectrum (F), we derive 
the relative flux at the Nat D lines by comparing fluxes inside 
passbands in the center, red and blue part of each lines (see 
Sect.|2.5|and 3.1 for the description of the passbands): 


Trelit, A/1) = 


2 X FiO 


FiB) + F{R) 


(5) 


We then normalized the relative time serie to unity. For the 
smaller passbands, we average, spectrum by spectrum, the rel¬ 
ative fluxes of the two Nat D lines. 

Because we always compare together same parts of each 
spectrum, the relative flux should be constant with time except 
if absorption by the planetary atmosphere is present during the 
transit. With this method, we cannot apply any planetary radial- 
velocity correction. Taking into account the radial-velocity effect 
would change the central passband, thus it would change the rel¬ 
ative flux (not because there is an absorption, but because we 
do not have a constant central passband). However this is poten¬ 
tially an issue only for the smallest passbands for which we can 
expect a loss of absorption. 

Another type of telluric correction is applied to the rela¬ 
tive flux sequence. As the telluric absorption changes with air- 
mass, it is expected to see correlations between relative fluxes 
and the airmass, especially when the passbands contain telluric 
lines. The first correction applied on individual stellar spectra 
(Sect. |2.3| ) should have already mitigated this effect. However, 
some light curve still show correlations with airmass. To remove 
any residual airmass effect, we model the flux variation as a lin¬ 
ear function of the airmass and remove the linear trend. If possi¬ 
ble, we only consider out-of-transit data to perform our fit. 

Astudillo-Defru & Rojo|(2013|l discussed whether or not flt- 


6(AA) = 


T'reABn) 

T^reli^Out} 


- 1 . 


( 6 ) 


2.7. Systematic effects 

To measure the uncertainty on the relative absorption depth 
d(A/l), we propagated the errors of the reduced spectra through 
our analysis. We estimated the errors on the reduced spectra as 
random photon noise obeying Poisson statistics. Systematic ef¬ 
fects, on the other hand, are likely to contribute to the total noise 
budget. We have used different statistical methods, such as em¬ 
pirical Monte-Carlo ( Redfleld et al.|[2008] l or bootstrapping, to 
estimate the impact of correlated noise. The basic principle of 
these methods is to randomize data, to artificially create new set 
of observations, and to feed them to our data analysis. 

Firstly, we followed the empirical Monte-Carlo technic pre¬ 
sented in Redfleld et al.|P008] l. The idea is to see whether or not 
the measurement is an artifact of the data or if it is really created 
by the transit and hence due to the planetary companion. We se¬ 
lect a sub-sample of spectra of all the spectra available in a night. 
The choice of this sub-sample is then fixed (we call it a scenario). 
Next, this new sample is randomly divided in two parts to create 
an in-transit and an out-of-transit simulated data sets. These two 
data sets always contain the same number of spectra. 


As in Redfleld et al. (2008 1 , we explored three scenarios per 
night. In the first scenario (out-out), the spectra are only selected 
among the nominal out-of-transit spectra. The amount of spec¬ 
tra in the simulated in-transit data compared to the number of 
spectra in the simulated out-of-transit data are in the same pro¬ 
portion as in the nominal observations. As we choose randomly 
only among the out-of-transit spectra, we expect to And the dis¬ 
tribution of computed relative absorption depth centered at zero. 
Similarly, an in-in scenario is created, but the random in-transit 
and out-of-transit data are chosen only in the spectrum observed 
during the transit. As the planetary signal should be present in¬ 
side every in-transit observation, the final in-in distribution is 
also expected to be centered at zero. Finally an in-out scenario is 
executed. All the out-of-transit observation form a fixed master- 
out (actually the same than the nominal master-out). We choose 
randomly a sub-sample of in-transit spectra. Here, we did not 
fix the number of spectra of the in-transit sub-sample. This later 
one is composed by half to the totality of the number of nominal 
in-transit spectra. 


How to interpret the computed distribution? [Redfleld et al. 
( |2008] l considered the standard deviation of the out-out distribu¬ 
tion as the global error on the transmitted signal. Nonetheless, 


ting a transit model to the data can change the measurement of 
the absorption depth. They showed that the absorption depth is 
well measured by simple average differences. Indeed, differen¬ 
tial stellar limb-darkening do not affect significantly measure¬ 
ment in such narrow passbands. The transmitted flux or relative 
depth is then given by the difference of the average of the relative 
flux in and out-of-transit: 


Astudillo-Defru & Rojo (20131 considered instead the standard 


deviation of the in-out distribution as the error on the transmitted 
signal. Here, we follow Redfleld et al.| ( |2008] l and take the out- 
out scenario to infer the error on our measured absorption depth, 
because it is independent of any planetary signal. Nevertheless, 
the errors are overestimated with the standard deviation of the 
out-out scenario, because in each iteration only a fraction of the 
data is used. Thus, to get correct errors on our absorption depth, 
we divide the standard deviation of the out-out scenario by the 
square root of the ratio of the total number of spectra to the num¬ 
ber of spectra out-of-transit. 


Integrating the transit light curve over the in-transit dura¬ 
tion makes it possible to compare the result with the binned 
atmospheric absorption depths described in Sect. [2.5 j Both ap¬ 
proaches give similar results, because of the normalization pro¬ 
cesses and the use of photon noise weighted averages in the two 
respective methods. 


Secondly, we calculate false alarm probabilities. We random¬ 
ize spectra and recreate transmission spectra and light curves 
(i.e. by randomizing the time sequence), and quantify how many 
randomized data sets yield a signal as significant as the real mea¬ 
surement. This is done night by night and by considering all the 
original spectra. 
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Fig. 2. HARPS transmission spectrum of HD 189733b around the sodium Na i D doublet. Upper panel; Master out spectrum (stack 
of the out-of-transit spectra) normalized to unity. A telluric spectrum is also shown with a vertical offset to be able to identify 
water telluric lines. Middle panel: Overall transmission spectrum of the exoplanet atmosphere taking into account all the observing 
nights (light grey). The transmission spectrum is shown binned by 20x with black circles. The planetary radial-velocity correction 
is applied. We can easily see the two Nai D lines core absorption from the planetary atmosphere. We show a gaussian fit to each 
Nai D lines (red). We measure line contrasts of 0.64 + 0.07% (D2) and 0.40 + 0.07% (Dl) and FWHMs of 0.52 + 0.08 A. The 
transition wavelengths of the doublet, in the planet rest frame, are indicated with blue dashed line. A net blueshift is measured, by 
our ht, with a value of 0.16 + 0.04 A. This blueshift may be due to wind in the upper layer of the atmosphere with a velocity of 8 + 2 
km/s (see text). Lower panel; Residuals to the gaussian fit. The noisy part is centered on the two sodium lines wavelength where the 
stellar flux is the lowest. 


3. Results and analysis 

3. L Na I D detection 

As we focus on the two transitions of the sodium doublet (D2 at 
45889.951 & Dl at 45895.924 A), we restrict the wavelength 
range for our analysis to the domain from 5870 to 5916 A. 
Taking into account the systemic velocity of -2.2765 km s~', 
these lines are, in the solar system barycentric rest frame, at 
45889.906 and 45895.879 A, respectively. 

To compute absorption depths, we define different cen¬ 
tral passbands (C) composed by two sub-bands centered 
on each sodium transition lines. Total band width of 
6=2x3 A, 3=2xl.5 A, 1.5=2x0.75 A, 0.75=2x0.375 A and 
0.375=2x0.188 A were chosen (the 6=2x3 A is actually 
2x2.98 A). We add to these bands a larger passband of 12 A 
to perform the comparison of our results to those of |Redheld| 


|et al.| ( |20081 l; [Jensen et al.| ( |201 Ij l; [Huitson et al.j ( |2012] l. As 
the 12 A band encompasses the two sodium lines, we adjusted 
it on the center of the doublet. For every bands, we choose a 
unique reference passband corresponding to the blue (B) and 
red (R) adjacent passbands of the central passband of 12 A (i.e. 
B = 5874.89 - 5886.89 A, R = 5898.89 - 5910.89 A). 

3.2. Transmission spectrum anaiysis 

The total signal-to-noise ratio (SNR) per extracted pixel of the 
in and out master spectra (7y„ and 'Tout) ranges from 400 to 550 
for all nights, in the continuum around 589 nm. A pixel on the 
HARPS detector represents ~ 0.8 km s ' (or 0.016 A at 589 nm). 
The total SNR for each master is about 850 when co-adding the 
nights. Thus, a hnal SNR on the transmission spectra (91) of 
about 1 200 per extracted pixel is reached. For the transmission 
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Fig. 3. HARPS transmission light curve of HD 189733b. The three nights are averaged together with bins of 3 spectra. Left; The 
relatives fluxes are integrated in each spectrum over a 12 A passband including the whole sodium Nat D doublet and compared to 
reference passbands in the same spectrum. The overall absorption depth is detected at 620 ppm with 7.8 cr. Right: The relatives 
fluxes are integrated over two passbands of 0.75 A around each sodium Nat D lines. The absorption depth is detected at 3 250 ppm 
with 9.8 cr. Note the different vertical scales. 


spectra, the SNR in the core of the sodium lines decreases to 
~ 260. It is therefore important to co-add the fluxes over tens to 
hundreds pixels to reach a sufficient SNR for atmospheric detec¬ 
tions. We compute different transmission spectra, notably with 
and without planetary radial-velocity correction. The standard 
deviation of the transmission spectra was about 2 200 ppm for 
each night in the continuum. When adding the nights together 
the standard deviation goes down to 1450 ppm. As we can see 
in Fig. 1^ the two exoplanetary sodium lines peak out from the 
continuum. Without planetary radial-velocity correction, the at¬ 
mospheric absorption is best detected in our 1.5=2x0.75 A pass- 


band at a level of 0.309+0.034% (9 cr). Our detection in the 12 A 
passband is measured at a level of 0.060 + 0.008% (7.5 cr). Here 
and in the following, we dehne absorption as a negative rela¬ 
tive depth. With planetary radial-velocity correction the signal in 
these same 1.5=2x0.75 and 12 A passbands are 0.320 ± 0.031% 
(10.3 cr) and 0.056 ± 0.007% (8 cr) respectively (see Table [^. 
As the radial-velocity correction consists of a wavelength shift 
smaller than 0.25 A, it is normal that the integrated signal on 
these passbands are consistent with each other. Furthermore, our 
measuremen ts are in agreement with the detection by [Huitson 


et al. 


_ ^(20121 from space (0.051 ± 0.006% for the 12 A pass- 

bands) and by [Redheld et al. ( 2008] l and Jensen et al. ( 201 l| l 


from ground (0.067 + 0.020 and 0.053 + 0.017 for the 12 A pass- 
band, respectively). 


These absorption signals can be interpreted as equivalent 
relative altitudes. Altitudes can be inferred by considering the 
transmission at the limb. Indeed, the atmospheric absorption 
(which depends on wavelength) is due to an optically thick layer 
(presence of an absorber) at a certain height above the mea¬ 
sured planetary radius (in broad band light curve). The height 
of this layer is generally a few times bigger than the atmospheric 
scale height of the planet {H = kTlug), because of the variation 


with wavelength of the cross section of the absorber (Lecavelier 
Des Etangs et al. 2008a|l. The absorption depths measured in 


the 1.5=2x0.75 A passband correspond to equivalent altitude 
of 5100 ± 500 km assuming unresolved features. This is about 
27 atmospheric scale height of HD 189733b {H = 190 km for 


1 140 K, the equilibrium temperature of HD 189733b, assuming 
a mean molecular weight of /i = 2.3 and an albedo of 0.2). 

In order to understand the impact of our radial-velocity cor¬ 
rection it is then worth to look at the smallest passband. At hrst 
glance, integrations over 0.75=2x0.375 A and 0.375=2x0.188 A 
passbands give us smaller signals than the signals we acquire 
without radial-velocity correction. This is not expected, since 
the radial-velocity correction should add up all the absorptions 
signals and thus strengthen the signal. A closer look to the trans¬ 
mission spectrum with a Gaussian htting to each line shows that 
the exoplanetary sodium lines are blueshifted by about 0.16 A 
in the planet rest frame. Integrating the lines including this shift 
(with the smallest passbands blueshifted) enables the full signal 
to be recovered. 


What can produce the blueshift of the exoplanetary sodium 
lines we measure and is it signihcant? A global error on the plan¬ 
etary radial-velocity shift is likely to be lower than 3 km s ' 
(0.06 A) considering the HARPS precision on the stellar radial 
velocity propagated to the planet. Our 0.16 A value is by fa r 

(12008 I. 


Redheld et al. 


smaller than the ~ 0.75 A measured by 
Since this specihc blueshift would have been easily measured 
with our data, we do not conhrm the value reported by these 
authors. Nevertheless, our detected 0.16 + 0.04 A shift is sig¬ 
nihcant and most likely real, since it is seen even on the trans¬ 
mission spectra without radial velocity correction. A possible 
physical explanation is that a net blueshift is imprinted by winds 
within the exoplanet atmosphere (see [Snellen et al.||20f0| for 
a hrst detection of winds in exoplanets). Theoretical models 
{e.g. [Kemp ton et al.j 2014| Showman et al. 2013| ) indeed pre¬ 
dict a net blueshift due to both planetary rotation and winds in 
the terminator of the atmosphere (the part we are probing with 
transmission spectroscopy). The blueshift we measured corre¬ 
sponds to a wind speed of 8 + 2 km/s. Models typically pre¬ 
dict a 3 km/s blueshift. However, these winds are estimated for 
pressures down to ~ 10“^ - 10“® bar. As we will see later, our 
sodium detection allows us to probe pressures of ~ 10“^ - 10“® 
bar (0.1 - 0.001//bar) at altitude of 10000 km. We can therefore 
expect stronger winds than 3 km/s in the higher atmosphere. This 
calls for additional theoretical works on winds in the upper at- 
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Fig. 4. Distributions of the empirical Monte-Carlo analysis for the 12 A passband. Top row; results for the transmission spectrum 
method. Bottom row; results for the light curve method. The columns correspond to the three different nights of observation (nights 
1 to 3, from left to right), in-in, out-out and in-out scenarios are shown in blue, red and green. Distributions are shown in continuous 
lines, Gaussian fit to the distribution in dotted lines and vertical lines are the average values. Results on the measurements and errors 
coming from these distributions are summarized in Fig.j^where we can compare it to the nominal values. 








Fig. 5. Distributions of the false-alarm probability analysis for the 12 A passband. Top row; results for the transmission spectrum 
method. Bottom row; results for the light curve method. The columns correspond to the three different nights of observation (night 
1 to 3, from left to right). Distributions are shown in continuous lines, gaussian fit in dotted lines and nominal values to infer 
false-alarm probabilities in vertical lines. 


mospheres of exoplanets, and new data to confirm the blueshift 3.3. Transmission iight curve anaiysis 
value. 

We compute the transmission light curve of HD 189733b for 
every observation nights following Sect. |2.6| and for every 
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Fig. 6. Comparison of the different methods. The relative absorption depth due to sodium lines are shown for all the nights and 
for the different methods. Grey squares are the measurement levels without any correction. Particularly, telluric lines imprint spu¬ 
rious systematics in any directions. Magenta circles and diamonds are measurements made with the light curve method and the 
empirical Monte-Carlo. Blue left triangle and diamond are for the transmission spectrum without radial-velocity correction and 
the corresponding empirical Monte-Carlo. Cyan right triangles and diamonds are for the transmission spectrum with planetary 
radial-velocity correction and the corresponding empirical Monte-Carl o. Total absorption de pth when adding all the night are also 


indicated. Left: Summary for the 12 A passband. The measurement of 


Huitson et al. 


(2012 1 (3 transits with HST/STIS) is shown 


through black hexagon. The measurement of Redfield et al.| ( |20()8) l (35 exposures in-transit (600 s) observed during 11 transit events) 
is shown through black left triangle and diamond. The errors correspond to their photon noise and empirical Monte-Carlo analysis, 
respectively. HARPS precision is comparable with HST/STIS precision showing the importance of ground-based high resolution 
observation of exoplanetary atmospheres. Right: Summary for the 1.5 -2x 0.75 A passband. 


Table 3. Summary of the measured relative absorption depth in [%] on the transmission spectra with radial-velocity correction. 


AA [A] 

0.375=2x0.188 

0.75=2x0.375 

1.5=2x0.75 

3=2xl.5 

6=2x3 

12 

# pixel 

~ 2x 12 

~ 2x 25 

-2x50 

-2x98 

- 2 X 194 

- 780 

Night 1 
Night 2 
Night 3 

0.764±0.095 

0.495±0.112 

0.560±0.104 

0.532±0.064 
0.431 ±0.077 
0.500±0.071 

0.378±0.046 

0.198±0.054 

0.415±0.050 

0.202±0.025 
0.028±0.029 
0.221 ±0.027 

0.115±0.014 

0.015±0.017 

0.113±0.015 

0.077±0.010 

0.018±0.011 

0.080±0.010 

All Nights 

0.571±0.065 

0.472±0.044 

0.320±0.031 

0.141±0.017 

0.075±0.010 

0.056±0.007 


Table 4. Summary of the measured relative absorption depth in [%] on the light curve. 


AA [A] 

0.375=2x0.188 

0.75=2x0.375 

1.5=2x0.75 

3=2xl.5 

6=2x3 

12 

# pixel 

- 2x 12 

- 2 X 25 

-2x50 

-2x98 

- 2 X 194 

-780 

Night 1 
Night 2 
Night 3 

0.606±0.179 
0.244±0.199 
0.521 ±0.224 

0.468±0.095 

0.165±0.106 

0.661±0.119 

0.339±0.053 

0.143±0.059 

0.457±0.066 

0.207±0.031 

0.040±0.034 

0.274±0.039 

0.127±0.018 

0.027±0.021 

0.133±0.023 

0.064±0.013 

0.026±0.014 

0.091±0.016 

All Nights 

0.478±0.113 

0.451 ±0.060 

0.325±0.033 

0.184±0.019 

0.102±0.012 

0.062±0.008 


passbands described in Sect. 3.1 For our analysis, the transit 
ephemeris are fixed. They are given by the parameters described 
in Table 1^ For the 12 A passband, we estimate errors on each 
individual spectrum of 250 ppm to 800 ppm. The propagation of 
these errors to the time series provides us uncertainties of about 
100-150 ppm on the values of the baseline and on the absorp¬ 
tion depth. This allows us to determine absorption depths due to 
sodium in the exoplanet atmosphere during each night at a level 
between 1.9 and 5.7 cr (see Table [^. Note that when we sub¬ 
tract the best-fit absorption depth, the standard deviation of the 
residuals is about 550 ppm for each night. When co-adding data 
of the three nights together, we obtain an absorption signal of 


0.062 ± 0.008% (7.8 cr) for the 12 A passband (Fig.[^. Our best 
measurement is achieved for the passband of 1.5 = 2 x 0.75 A 
and delivers a value of 0.325 + 0.033% (9.8 cr). This is in perfect 
agreement with our transmission spectrum method and shows 
both the complementarity and the compatibility of the two meth¬ 
ods. Compatibility because the measured absorption level are 
within a 1 cr error bars (see Astudillo-Defru & Rojo|20T3 i and 
complementarity because the second method carries line infor¬ 
mation and demonstrates that the additional absorption is indeed 
observed during the transit (see also Sect. |3.4| l. As expected, a 
small decrease on the measured absorption signal is observed 
for the smallest passband (0.375=2x0.188 A), which width is 
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commensurable with the wavelength shift imprinted by the ra¬ 
dial velocity of the planet. This effect thus results in a flux loss, 
hence to a weaker signal. 


3.4. Systematic effects 


Possible systematic errors can be present in the data. To esti¬ 
mate such effects, we perform an empirical Monte-Carlo (EMC) 


analysis for each of our observed transits following Sect. 2.7 


This analysis makes us confident that our measured signals is in¬ 
deed due to the exoplanetary transit and not spurious (e.g. due 
to other astrophysical effects as stellar rotation, observational 
conditions, or instrumental effects). The EMC is performed for 
both the transmission spectrum and the light curve methods us¬ 
ing 3 000 iterations for each scenario described in Sect. |2.7| (in- 
in,out-out,in-out) and for each night. The standard deviation of 
the distributions does not increase significantly when increasing 
the number of iterations. 

The results are shown in Eig. and the errors are summa¬ 
rized with all the different methods in Eig. We can see that 
all the in-in and out-out distribution are centered on zero show¬ 
ing that the measured absorption signal comes indeed from the 
transit and has no other sources of explanation. The in-out dis¬ 
tributions are well centered at the value determined by the trans¬ 
mission spectrum and light curve methods. The errors given by 
the EMC are larger than the nominal errors (yielding from the 
propagation of the photon noise), but only by a factor of < 2. 
This shows that while systematic errors are present in our data 
they do not undermine our detection. The empirical Monte-Carlo 
simulation described above strengthens the transit origin of the 
detected signal. 

We also computed false alarm probabilities for each night for 
both the transmission spectrum and the light curve methods. We 
performed 10000 iterations and found false-alarm probabilities 
to our sodium detection of ~ 0.5%, ~ 9% and < 0.01% for night 
1, 2, and 3 respectively (see Eig. Adding all the nights yields 
a false-alarm probabilities < 0.01%. 


4. Discussion 

4.1. Summary of the Na i D detection 

A summary of our results is shown in Eig. All the different 
methods to measure the transit depth i5(A/l) in a given wave¬ 
length bin AT are shown for each night individually and also 
for the combined data set, for the 1.5 = 2 x 0.75 A and 12 A 
passbands. 

A number of aspects must be pointed out: Eirst, the correc¬ 
tion of the telluric lines is mandatory. In fact, water lines around 
the sodium Nai D doublet imprint spurious systematics in any 
possible direction. It depends if the transit is observed at low or 
high airmass compared to the out-of-transit data and at which 
value of the barycentric Earth radial velocity (BERV). Indeed, 
due the BERV, the telluric lines will not always fall at the same 
place in the stellar spectra. Eor a fixed choice of passbands, a 
given telluric line may not be always included. Then, the ratio 
with the reference passbands, which contain other telluric lines, 
can be affected in any direction. Eor example, during night 3, 
the signal measured without telluric correction is, by chance, the 
same signal measured after telluric correction, as the positions 
of telluric lines in the transmission spectrum do not impact the 
sodium signal extraction. It should be noted that the correction 
for night 2 is less robust than the one of the two other nights. 



Fig. 7. Comparison of the transmission spectra around the Nat 


doublet obtained by the HST/STIS instrument (Huitson et al. 
2012[ in blue) and by the HARPS spectrograph (binned by 5x, 


in black). 


because of the lack of a baseline before the transit. This im¬ 
plies that we probably over-corrected the data taken during this 
night, hence the measured value is under estimated. Secondly, 
the absorption depths measured with transmission-spectra and 
light-curve methods are consistent with each other within the 
same night. Thirdly, systematic effects are different from night 
to night. Nights 2 and 3 were probably more stable than night 
1. Systematics in night 1 are probably also bigger due to the 
low-cadence of the observations. Even if similar signal-to-noise 
ratios are reached in the three different nights, it seems that it 
is better to have a high-cadence observational strategy (provided 
all the spectra have a SNR above a certain threshold). We also 
note that the measured values are different for different nights. 
At this level of precision, and assuming we are efficiently cor¬ 
recting for systematics, we cannot exclude intrinsic variability 
related to the stellar or planetary properties. Our results favor a 
variation from night to night within 1-2 cr. We did not take into 
account the potential effect of a differential limb darkening or 
of the Rossiter-McLaughlin effect. However, the impact of dif¬ 
ferential limb-darkening should be an order of magnitude lower 
than the achieved precision (Charbonneau et al.|2002[ [Redfield 


et al.|200^. On the other hand, the Rossiter-McLaughlin effect 


can mimic an exoplanetary atmospheric absorption and imprint 
features at similar level in a transmission spectrum. Nonetheless 
it should average out during a whole transit if the exposures are 
equally distributed over the transit ( |Dravins et al.|2014) l. This is 
the case for our three analyzed transit events. Einally, we would 
like to point out that the precision obtained from the coaddi¬ 
tion of all the nights is comparable with HST/STIS precision 
(with the same number of observed transits) or with 10-m sizes 
telescopes. This highlights the potential of ground-based high 
resolution observations of exoplanetary atmospheres from 4-m 
telescopes. 

As mentioned before, we obtain the same precision on our 
detection (0.056 + 0.007%) as on the HST/STIS measurements 
of Huitson et al. ( 20I2| l (0.051 + 0.006%) on a relatively wide 
passband (here the 12 A). We can take advantage of the HARPS 
high spectral resolution (7? ~ 115 000), compared to HST/STIS 
(7? ~ 5 500) to characterize the narrow line cores of the Nat 
doublet and, as a consequence, exploring higher regions in the 
atmosphere of HD 189733b. 
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Fig. 8. Fit of a range of tj models to the HARPS transmission spectra of HD 189733b. The data and the models are binned by 20x 
(panel a), 5x (panel b) and 2x (panels c and d). Models correspond to those described in Table|^ Some additional models (2100 K, 
2600 K, 3100 K, the hottest in panels a and b) and the continuum level (in red), which correspond to the altitudes of a haze layer 
deck, are also shown for comparison. Our analysis shows that the atmosphere of HD 189733b is non-isothermal. 


Table 5. Choice of the wavelength ranges for the atmospheric model fits, number of Nai lines FWHM encompassed in the wave¬ 
length ranges, corresponding altitudes, and fitted temperatures. We centered the wavelength ranges at the center of each sodium line 
in the planet rest frame taking into account the observed blueshift. 


Model Atmospheres 

Wavelength Ranges [A] 

#NaiD’sFWHM 

Corresponding Altitude [km] 

Fitted Temperature [K] 

Model 0 

5870.00-5882.22 

5903.24-5916.00 


0 

1 140 

Model 1 

(Wing' shoulders) 

5882.22-5889.22 

5890.26-5895.20 

5896.24-5903.24 

2-14 

1500±1500 

1630±70 

Model 2a (D1 core) 

5895.20-5895.46 

5895.98-5896.24 

1-2 

2700±800 

1700±320 

Model 2b (D2 core) 

5889.22-5889.48 

5890.00-5900.26 

1-2 

3 800±900 

2170±320 

Model 3a (D1 core) 

5895.46-5895.67 

5895.77-5895.98 

0.2-1 

5100±3 100 

2220±340 

Model 3b (D2 core) 

5889.48-5889.69 

5889.79-5890.00 

0.2-1 

7900±5500 

3 220±270 

Model 4a (D1 core) 

5895.67-5895.77 

^0.2 

9800±2800 

2600±600 

Model 4b (D2 core) 

5889.69-5889.79 

<0.2 

12700±2600 

3 270±330 


The presence of haze in the HD 189733b atmosphere has 
been previously deduced from low-resolution spectra covering 
much broader spectral regions than the one considered here (e.g. 
Lecavelier Des Etangs et al. 2008a Pont et al. 2013)1. Rayleigh 


in 


scattering by hazes dominates the transmission spectrum in the 
optical, only allowing the narrow core of the sodium lines to 
be observed at high spectral resolution. This is indeed what we 
observe in our data. 


As described in Sect. 3.1 we compute different absorption 
depths (5(A/l) for different passbands (summary in Table[^and|^. 
The absorption depth increases for narrow cores passbands. This 
shows that most of the signal we measure originates in the nar¬ 
row line cores (< lA), even if a small portion of the absorption 
comes from a broader component. Actually, a Gaussian fitted 
to each lines yields a full width at half maximum (FWHM) of 
0.52 + 0.08 A (see Fig. |^. A comparison to the resolution ele- 
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ment of 0.05 A (2.7 km ) shows that the lines are resolved by a 
factor of ~ 10 (see Fig0. The same fit allows us to measure ab¬ 
sorption depths in the core of the two Na i D lines of 0.64+0.07% 
and 0.40 ± 0.07% for the D2 and D1 lines, respectively. 


4.2. The Na i D sodium doublet as a probe to the upper 
atmosphere 


We compare our measurements with transmission spectroscopy 
models. We use the 77 model described in |Ehrenreich et al.] ( |2006| 
|2012| |2014[ ). The model computes the opacity along the line 
of sight grazing the atmospheric limb of planet and integrate 
over the whole limb. The planetary atmosphere is modeled as 
a perfect gas in hydrostatic equilibrium composed by 93% of 
molecular hydrogen and, 7% of helium, and solar abundance of 
atomic sodium (volume mixing ratio of 10“®). Several model at¬ 
mospheres are simulated with isothermal temperature profiles 
for temperature ranging from 1 000 K to 3 600 K by step of 
100 K. To increase the resolution in temperature, we interpo¬ 
late linearly between the models. The sodium line profiles are 
modeled with Voigt functions resulting from the convolution 
of a Doppler thermal profile (dominating the line cores) and 
a Lorentzian profile accounting for the natural and collisional 
broadenings in the line wings. Collisional or pressure broaden¬ 
ing is completely dominating the line shapes far (over scales of 
~ 100 nm) from the cores ( Iro et all|20051 l. The half-width at 
half-maximum (HWHM) of the pressure broadened Lorentzian 
follows the prescriptions of Burrows et al.| ( 2000)l and Iro et al. 


( 2005| l, which are valid over the wavelength range we are study¬ 
ing here. The modeled transmission spectra, calculated at a res¬ 
olution of O.OIA, are convolved with the average HARPS Line 
Spread Lunction (LSL), which is well represented by a Gaussian 


with width of 0.05A. In contrast with HST/STIS, the comparison 
between the models and our data is eased by the fact that HARPS 
resolves individual lines arising in the planet atmosphere by a 
factor ~ 10 (see above). 


Before comparing our transmission spectrum to models, we 
first fitted to our data a simple haze model. Knowing that the 
Rayleigh slope in our wavelength range analysis is negligible, 
hazes can be well represented by an absorption cut-off at a given 
level. This level must be compatible with our measured contin¬ 
uum, which corresponds to the wavelength ranges of model 0 
in Table We want also to determine if the sodium wings ab¬ 
sorption dominates the haze absorption in the transmission spec¬ 
trum. The wings’ shoulders region is described by the wave¬ 
length ranges of model 1 in Table Therefore, the maximum 
altitudes of the cloud deck which is compatible with the contin¬ 
uum level is given by a constant fit to the wavelength regions 
used to fit model atmospheres 0 and 1 to the data. This constant 
level gives a Bayesian information criterion (BIC) of 1 597 in the 
wings region (region 1 ). 

To compare our observed transmission spectrum to models, 
our approach is to fit the different regions of the measured trans¬ 
mission spectrum around the sodium doublet with model atmo¬ 
spheres with different isothermal profiles. We divided the trans¬ 
mission spectrum in separated wavelength ranges around each 
line of the sodium doublet, each corresponding to altitude slices. 
The altitude slices are given by the maximum between the alti¬ 
tude errors on the data and the difference in altitudes in the given 
ranges. These ranges are listed in Table We fit these differ¬ 
ent parts of the spectrum with different isothermal models, prob¬ 
ing separately the line cores (models 4), the region encompassed 
within IxLWHM (excluding the line cores; models 3), the re¬ 



Fig. 9. Temperature profile obtained by our fitting process of rj 
models to the transmission spectrum. The right and left triangle 
correspond to the Nai D1 and D2 lines, respectively. The dia¬ 
mond correspond to the wing’shoulders. A temperature gradient 
of ~ 0.2 K km“' is measured. 


gion between 1 x and 2x the LWHM (models 2), the line wings 
(excluding the previous regions; model 1 ), and the continuum 
(model 0). All the fits are obtained with a^^ minimization over 
the grids of atmospheric models obtained by varying the tem¬ 
perature and a general offset in relative flux. In the continuum 
wavelength range we cannot see a significant sodium absorp¬ 
tion due to absorbing wings. Then, model 0 was chosen with 
a fixed temperature of 1140 K (the equilibrium temperature of 
HD 189733b). The fit of this model to the continuum allows us to 
determine the offset in planetary radius between the data and the 
models. We fix this offset for the subsequent model adjustments. 
We fit models 1, 2, 3, and 4 to their respective spectral regions, 
simply considering the model temperature as a free parameter. 
The fit of model 1 gives us a temperature of 1 630 K for a BIC 
of 1 542. The comparison with the BIC of 1 597 obtained with a 
constant haze level is very strong evidence (ABIC > 20) in favor 
of a scenario where the Na i wings absorption dominates the haze 
absorption. Thus, part of our measured sodium absorption is due 
to the wings and correspond to lower altitudes and temperatures. 
The three other models (2, 3 and 4) investigate different parts 
of the line cores inside 1 A (2 LWHM, see Table |^. The best 
fit models are shown in Lig. Lor each range of altitudes, we 
therefore derive one temperature. The resulting temperature pro¬ 
file as a function of altitude is shown in Lig. The temperature 
linearly increases with altitude with a gradient of ~ 0.2 K km ' 
(0.2 ± 0.1 K km“'). Note that our measured temperatures are un¬ 
derestimated due to the use of isothermal models, which give, 
for a same upper temperature, more extended atmospheres than 
models with positive temperature gradient. 

Our high-resolution measurement of sodium in HD 189733b, 
allows us to probe a new atmospheric region (between ~ 1.1 - 1.2 
planetary radius, corresponding to pressures of 10 “^ - 10 ® bar), 
above the result previously obtained at lower resolution ( |Huitson| 
et al.|2012] l. The present study suggests that a part of the sodium 
absorption takes place up to the (lower) thermosphere, where 


heating by the stellar X/LUV photons occurs (Lammer et al. 

2003 

Yelle 2004 |Vidal-Madjar et al.|2011| Koskinen et al. 

2013 

1 . Meanwhile our interpretation relies on hydrostatic mod- 


els, which could be questionable at very high altitudes. Lurther 
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theoretical analysis of these results is developed in a separate 
paper (|Heng et al.|2015[ accepted in ApJ). 


5. Conclusions 


We have presented an analysis of all the available transit obser¬ 
vations of HD 189733b with HARPS on ESO 3.6 m telescope. 
We carefully corrected for the change in radial velocity of the 
planet and for telluric contamination, fully exploiting the high 
spectral resolution and the stability of HARPS. We detect excess 
of absorption due to Nai D lines in our transmission spectrum. 
We measure line contrasts of 0.64+0.07% (D2) and 0.40+0.07% 
(Dl) and FWHMs of 0.52 + 0.08 A. Sodium is clearly detected 
over several different passbands and the signatures are robust 
against different statistical tests. Moreover, these detections are 
consistent and comparable in terms of precision with those ob¬ 
tained from space-borne or 10 m class ground-based facilities. 

High-resolution permits us to measure a blueshift in the line 


positions of 0.16 + 0.04 A. We interpret it as wind in the upper 
layer of the atmosphere of a velocity of 8 + 2 km/s. Since we 
resolved the sodium lines, we used it as a probe of the lower 
pressure part of the atmosphere. We measured a non-isothermal 
temperature profile inside our range of probed altitudes, where 
the sodium absorption is peaking out above the haze present in 
the atmosphere. We found that the temperature increases towards 
the thermosphere. A temperature gradient of ~ 0.2 K km ' is 
measured. The existence of a positive temperature gradient does 
not depend on assumptions about the haze layer. In fact, the gra¬ 
dient can be measured from the resolved line cores and in partic¬ 
ular from the difference in absorption levels between the Dl and 
the D2 lines. This shows that beyond the detections of atomic 
species, transit spectroscopy allows characterization of physical 
conditions present in atmospheres. 


et al. 


These transit spectroscopy data were already used by Triaud 
( 2009| l and by [Collier Cameron et al. ( 2010) 1 to study the 


Rossiter-McLaughlin effect. Furthermore, with these same ob¬ 
servations, |DrGlo^(|20^ to be submitted) detected a slope in 
the planet-star radius ratios that Pont et al. (2008 20131 inter¬ 
preted as Rayleigh scattering. These studies demonstrate the rel¬ 
evance of studying exoplanet atmospheres with high-resolution 
spectrographs mounted on 4-meter-class telescopes. This is es¬ 
pecially important considering that several such facilities will 
be built in the coming years to prepare for the follow-up of ex¬ 
oplanet candidates from incoming space missions such as K2 
(on-going), TESS, CHEOPS, and PLATO. 

In the near future, the HARPS-like instrument ESPRESSO, 
will be mounted on the Very-Large Telescope (VET) incoher¬ 
ent focus. ESPRESSO could be used to probe exoplanetary at¬ 
mospheres in the optical. This will bring a more complete cen¬ 
sus of the chemical composition and study of aeronomic prop¬ 
erties of a large variety of exoplanets, including those with 
larger bulk densities or with fainter host stars. On the infrared 
side, several spectrographs will be built soon such as SPIROU, 
CARMENES, IRD and CRIRES+. Based on recent results {e.g. 
[Birkby et al.) 2013) [Snellen et al.||201^ a strong development 
of atmospheric detections and characterization in the infrared 
is foreseeable. Eventually, high-resolution spectra of exoplan¬ 
ets with broad spectral coverage will allow inventory of the at¬ 
mospheric chemical composition, the exploration of aeronomic 
properties such as pressure-temperature profile and lead to a bet¬ 
ter understanding of the formation and evolution of exoplanets. 
This is also a motivation for building high-resolution spectro¬ 
graph on the European Extremely-Large-Telescope e.g. HiReS 
and METIS ( [Udry et al.|2^ [Brandi et al.|20T4| [Snellen et al.j 


|2015[ ). In the meantime, extensive studies have to be done with 
existing facilities and data especially in the optical domain {e.g. 
Hoeijmakers et al.|2014 1 to help observers make the best use of 


future facilities. The present study is a signihcant step in this di¬ 
rection. It could be extended to the search of other species or as 
a systematic search in other hot-jupiter atmospheres. 
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